home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / lmfit.pro < prev    next >
Text File  |  1997-07-08  |  11KB  |  286 lines
  1. ; $Id: lmfit.pro,v 1.7 1997/02/13 01:57:43 stevep Exp $
  2. ;
  3. ; Copyright (c) 1988-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5.  
  6. function lmfit, x, y, a, weights=weights, fita=fita, $
  7.                     Function_Name = Function_Name, $
  8.                         alpha=alpha,covar=covar,$
  9.                         itmax=itmax,iter=iter,tol=tol,chisq=chisq, $
  10.                         itmin=itmin,double=double,$
  11.                         SIGMA=SIGMA,CONVERGENCE=CONVERGENCE
  12. ;+
  13. ; NAME:
  14. ;       LMFIT
  15. ;
  16. ; PURPOSE:
  17. ;       Non-linear least squares fit to a function of an arbitrary 
  18. ;       number of parameters.  The function may be any non-linear 
  19. ;       function.  If available, partial derivatives can be calculated by 
  20. ;       the user function, else this routine will estimate partial derivatives
  21. ;       with a forward difference approximation.
  22. ;
  23. ; CATEGORY:
  24. ;       E2 - Curve and Surface Fitting.
  25. ;
  26. ; CALLING SEQUENCE:
  27. ;       Result = LMFIT(X, Y, A, FITA=FITA, FUNCTION_NAME = name, ITER=ITER,$
  28. ;                ITMAX=ITMAX, SIGMA=SIGMA, TOL=TOL, YFIT=YFIT, WEIGHTS=WEIGHTS)
  29. ;                CONVERGENCE=CONVERGENCE
  30. ;
  31. ; INPUTS:
  32. ;   X:  A row vector of independent variables.  This routine does
  33. ;       not manipulate or use values in X, it simply passes X
  34. ;       to the user-written function.
  35. ;
  36. ;   Y:  A row vector containing the dependent variable.
  37. ;
  38. ;   A:  A vector that contains the initial estimate for each parameter.  
  39. ;
  40. ; WEIGHTS: 
  41. ;   Set this keyword equal to a vector of fitting weights for Y(i). This 
  42. ;   vector must be the same length as X and Y.  If instrumental (Gaussian)
  43. ;   weighting is desired, that is the measurement errors or standard deviations 
  44. ;   of Y are known (SIGMA), then WEIGHTS should be  set to 1/(SIGMA^2.) 
  45. ;   Or if statistical (Poisson) weighting is appropriate, set WEIGHTS=1/Y.
  46. ;   If WEIGHTS is not specified then, no weighting is assumed (WEIGHTS=1.0).
  47. ;
  48. ; FITA:
  49. ;   A vector, with as many elements as A, which contains a Zero for
  50. ;   each fixed parameter, and a non-zero value for elements of A to 
  51. ;   fit. If not supplied, all parameters are taken to be non-fixed.
  52. ;
  53. ; KEYWORDS:
  54. ;       FUNCTION_NAME:  The name of the function (actually, a procedure) to 
  55. ;       fit.  If omitted, "LMFUNCT" is used. The procedure must be written as
  56. ;       described under RESTRICTIONS, below.
  57. ;
  58. ;       ALPHA:  The value of the Curvature matrix upon exit.
  59. ;       CHISQ:  The value of chi-squared on exit
  60. ;       CONVERGENCE: Returns 1 if the fit converges, 0 if it does
  61. ;               not meet the convergence criteria in ITMAX iterations,
  62. ;               or -1 if a singular matrix is encountered.
  63. ;    COVAR:  The value of the Covariance matrix upon exit.
  64. ;       DOUBLE: Set this keyword to force the computations to be performed 
  65. ;               in double precision.
  66. ;       ITMAX:  Maximum number of iterations. Default = 20.
  67. ;       ITER:   The actual number of iterations which were performed
  68. ;       SIGMA:  A vector of standard deviations for the returned parameters.
  69. ;       TOL:    The convergence tolerance. The routine returns when the
  70. ;               relative decrease in chi-squared is less than TOL in an 
  71. ;               interation. Default = 1.e-6.
  72. ;       YFIT:   The fitted function evaluated at the input X values.
  73. ;
  74. ; OUTPUTS:
  75. ;       Returns a vector of calculated parameters, which are improvements
  76. ;       upon the initial guesses supplied in A.
  77. ;
  78. ; SIDE EFFECTS:
  79. ;       None.
  80. ;
  81. ; RESTRICTIONS:
  82. ;       The function to be fit must be defined and called LMFUNCT,
  83. ;       unless the FUNCTION_NAME keyword is supplied.  This function,
  84. ;       must accept a single value of X (the independent variable), and A 
  85. ;       (the fitted function's  parameter values), and return  an 
  86. ;       array whose first (zeroth) element is the evalutated function
  87. ;       value, and next n_elements(A) elements are the partial derivatives
  88. ;       with respect to each parameter in A.
  89. ;
  90. ;       If X is passed in as a double, the returned vector MUST be of
  91. ;       type double as well. Likewise, if X is a float, the returned
  92. ;       vector must also be of type float.
  93. ;
  94. ;       For example, here is the default LMFUNCT in the IDL User's Libaray.
  95. ;       which is called as : out_array = LMFUNCT( X, A )
  96. ;
  97. ;
  98. ;    function lmfunct,x,a
  99. ;
  100. ;         ;Return a vector appropriate for LMFIT
  101. ;         ;
  102. ;         ;The function being fit is of the following form:
  103. ;         ;  F(x) = A(0) * exp( A(1) * X) + A(2) = bx+A(2)
  104. ;         ;
  105. ;         ;dF/dA(0) is dF(x)/dA(0) = exp(A(1)*X)
  106. ;         ;dF/dA(1) is dF(x)/dA(1) = A(0)*X*exp(A(1)*X) = bx * X
  107. ;         ;dF/dA(2) is dF(x)/dA(2) = 1.0
  108. ;         ;
  109. ;         ;return,[[F(x)],[dF/dA(0)],[dF/dA(1)],[dF/dA(2)]]
  110. ;         ;
  111. ;         ;Note: returning the required function in this manner
  112. ;         ;    ensures that if X is double the returned vector
  113. ;         ;    is also of type double. Other methods, such as
  114. ;         ;    evaluating size(x) are also valid.
  115. ;
  116. ;        bx=A(0)*exp(A(1)*X)
  117. ;        return,[ [bx+A(2)], [exp(A(1)*X)], [bx*X], [1.0] ]
  118. ;    end
  119. ;        
  120. ;
  121. ; PROCEDURE:
  122. ;       Based upon "MRQMIN", least squares fit to a non-linear
  123. ;       function, pages 683-688, Numerical Recipies in C, 2nd Edition,
  124. ;    Press, Teukolsky, Vettering, and Flannery, 1992.
  125. ;
  126. ;       "This method is the Gradient-expansion algorithm which
  127. ;       combines the best features of the gradient search with
  128. ;       the method of linearizing the fitting function."
  129. ;
  130. ;       Iterations are performed until three consequtive iterations fail
  131. ;       to chang the chi square changes by greater than TOL, or until
  132. ;       ITMAX, but at least ITMIN,  iterations have been  performed.
  133. ;
  134. ;       The initial guess of the parameter values should be
  135. ;       as close to the actual values as possible or the solution
  136. ;       may not converge.
  137. ;
  138. ;       The function may fail to converge, or it can encounter
  139. ;       a singular matrix. If this happens, the routine will fail
  140. ;       with the Numerical Recipes error message:
  141. ;      
  142. ;
  143. ; EXAMPLE:  
  144. ;        Fit a function of the form:
  145. ;            f(x)=a(0) * exp(a(1)*x) + a(2) + a(3) * sin(x)
  146. ;
  147. ;  Define a lmfit return function:
  148. ;
  149. ;  function myfunct,x,a
  150. ;
  151. ;       ;Return a vector appropriate for LMFIT
  152. ;
  153. ;       ;The function being fit is of the following form:
  154. ;       ;  F(x) = A(0) * exp( A(1) * X) + A(2) + A(3) * sin(x)
  155. ;
  156. ;
  157. ;       ; dF(x)/dA(0) = exp(A(1)*X)
  158. ;       ; dF(x)/dA(1) = A(0)*X*exp(A(1)*X) = bx * X
  159. ;       ; dF(x)/dA(2) = 1.0
  160. ;       ; dF(x)/dA(3) = sin(x)
  161. ;
  162. ;        bx=A(0)*exp(A(1)*X)
  163. ;        return,[[bx+A(2)+A(3)*sin(x)],[exp(A(1)*X)],[bx*X],[1.0],[sin(x)]]
  164. ;     end
  165. ;
  166. ;   pro run_lmfunct
  167. ;         x=findgen(40)/20.        ;Define indep & dep variables.
  168. ;         y=8.8 * exp( -9.9 * X) + 11.11 + 4.9 * sin(x)
  169. ;         sig=0.05 * y
  170. ;         a=[10.0,-7.0,9.0,4.0]        ;Initial guess
  171. ;         fita=[1,1,1,1]
  172. ;         ploterr,x,y,sig
  173. ;         yfit=lmfit(x,y,a,WEIGHTS=(1/sig^2.),FITA=FITA,$
  174. ;                  SIGMA=SIGMA,FUNCTION_NAME='myfunct')
  175. ;         oplot,x,yfit
  176. ;         for i=0,3 do print,i,a(i),format='("A (",i1,")= ",F6.2)'
  177. ;  end
  178. ;
  179. ; MODIFICATION HISTORY:
  180. ;       Written, SVP, RSI, June 1996.
  181. ;           
  182. ;-
  183.        inexcept=!except
  184.        !except=1
  185.        eps=1e-7
  186.        on_error,2              ;Return to caller if error
  187. ;
  188. ;      Enable Math error trapping
  189. ;
  190.         j=check_math(0,1)
  191.  
  192.        ndata = n_elements(x)    ; # of data points
  193.        ma = n_elements(a)       ; # of parameters
  194.        if ma le 0 then begin
  195.            message, 'A must have at least ONE parameter.'
  196.        endif
  197.        nfree = n_elements(y) - ma ; Degrees of freedom
  198.        if nfree le 0 then message, 'LMFIT - not enough data points.'
  199. ;
  200. ;      Process the keywords
  201. ; 
  202.        if n_elements(function_name) le 0 then function_name = "LMFUNCT"
  203.        if n_elements(tol) eq 0 then tol = 1.e-9 ;Convergence tolerance
  204.        if n_elements(itmin) eq 0 then itmin= 5 ;Minimum # iterations
  205.        if n_elements(itmax) eq 0 then itmax= 50    ;Maximum # iterations
  206.        if (itmin ge itmax) then itmax=itmin
  207.        do_double = keyword_set(double)
  208. ;
  209. ;      Prepare the FITA vector
  210. ;
  211.        if n_elements(FITA) eq 0 then FITA = replicate(1, ma)
  212.        if n_elements(FITA) ne ma then $
  213.          message, 'The number of elements in FITA must equal those of A'
  214.        FITA=fix(FITA)
  215. ;
  216. ;      Prepare the SIG vector, this is sqrt(1/WEIGHTS)
  217. ;
  218.        if n_elements(WEIGHTS) eq 0 then WEIGHTS = replicate(1.,ndata)
  219.        if n_elements(WEIGHTS) ne ndata then $
  220.          message,'The number of elements in WEIGHTS must equal those of X and Y'
  221.        weights = weights > eps
  222.        SIG=1.0/WEIGHTS
  223. ;
  224. ;       If x or y or sig or a is double precision, set do_double to true
  225. ;
  226.        If (reverse(size(a))) [1] eq 5 or $
  227.          (reverse(size(x))) [1] eq 5 or $
  228.          (reverse(size(y))) [1] eq 5 or $
  229.          (reverse(size(sig))) [1] eq 5  then do_double=1
  230. ;
  231. ;
  232.        
  233.        if do_double then begin
  234.            chisq=0.0D
  235.            xx=double(x) & yy=double(y) & a=double(a)
  236.            ssig=abs(double(sig)) > 1E-12
  237.        endif else begin
  238.            chisq=0.0
  239.            xx=float(x) & yy=float(y) & a=float(a)
  240.            ssig=abs(float(sig)) > 1E-6
  241.        endelse
  242.  
  243. ;
  244. ;       Warning! The following call is to the actual NR recipies code.
  245. ;       Direct calls to this by the user are not supported, as this
  246. ;       function call will be removed in a future version of IDL
  247. ;
  248.        MRQMIN, function_name, xx, yy, ssig, ndata, a, fita, ma, covar, $
  249.          alpha, chisq, alambda, $
  250.          DOUBLE=do_double, itmin=itmin, itmax=itmax, $
  251.          tolerance=tol, niter=iter
  252. ;
  253.        convergence=1
  254.        if alambda lt 0 then begin
  255.            convergence=-1
  256.            message, 'Warning: Singular Covariance Matrix Encountered, LMFIT aborted.', /INFORMATIONAL
  257.        endif else begin
  258.            if (iter ge itmax) then begin
  259.                convergence=0
  260.                message, 'Warning: Failed to Converge.', /INFORMATIONAL
  261.            endif
  262.        endelse
  263. ;
  264.        diag=lindgen(ma)*(ma+1)
  265.        SIGMA=sqrt(abs(covar[diag]))
  266.        
  267.        if do_double then yfit=dblarr(ndata) else yfit=fltarr(ndata)
  268.        for i=0,ndata-1 do $
  269.          yfit[i] = (call_function(function_name,xx[i],a))[0]
  270.  
  271. ;
  272. ;    Check for Math Errors
  273. ;
  274.        errs=['Divide by 0','Underflow','Overflow','Illegal Operand']
  275.        j=check_math()
  276.        for i=4,7 do if ISHFT(j,-i) and 1 then $
  277.          message, 'Warning: '+errs[i-4]+ ' Occured.',/INFORMATIONAL
  278.        j=check_math(0,0)
  279. ;
  280. ;  return the result
  281. ;
  282.        !except=inexcept
  283.        return,yfit 
  284. ;
  285. END
  286.